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ABSTRACT 

In an attempt to determine the outer scale of turbulence driven by localized sources, such 
as supernova explosions in the interstellar medium, we consider a forcing function given by 
the gradient of gaussian profiles localized at random positions. Different coherence times of 
the forcing function are considered. In order to isolate the effects specific to the nature of the 
forcing function we consider the case of a polytropic equation of state and restrict ourselves to 
forcing amplitudes such that the flow remains subsonic. When the coherence time is short, the 
outer scale agrees with the half-width of the gaussian. Longer coherence times can cause extra 
power at large scales, but this would not yield power law behavior at scales larger than that 
of the expansion waves. At scales smaller than the scale of the expansion waves the spectrum 
is close to power law with a spectral exponent of —2. The resulting flow is virtually free of 
vorticity. Viscous driving of vorticity turns out to be weak and self-amplification through the 
nonlinear term is found to be insignificant. No evidence for small scale dynamo action is 
found in cases where the magnetic induction equation is solved simultaneously with the other 
equations. 
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1 INTRODUCTION 

Turbulence plays an important role in many branches of as- 
trophysics. During the past decade simulations have been em- 
ployed to address fundamental questions in turbulent stellar con- 
vection zones ( Spruit. Nordl und, & T itle 1990), acc retion discs 
( Balbus & Hawlev 1998), and interstellar t urbulence iKorpi et alJ 
1 1999a : .Balsara et alj.2004l : iDib et aljEoO^ . In these three cases 
the nature of energy injection is quite different; convective instabil- 
ity, magneto-rotational instability, and explicit driving through su- 
pernova explosions, respectively. Nevertheless, important insights 
into the nature of astrophysical turbulence have been obtained 
by studying turbulence forced on large scales by adding random 
time-dependent large scale perturbations to the velocity. This al- 
lows the study of the turbulent cascade of kinetic energy from 
large scales to small er scales, and eventually down to th e dissipa- 
tive scale (see, e.g.. lBrandenburg & Subramanianll2005l for a re- 
view). While the scale of the aforementioned instabilities can in- 
deed be large (comparable with the system size), this is not so 
evident in the case of random forcing by supernova explosions. 
The energy released by supernova explosions, and to a lesser ex- 
tend also by st ellar winds and outflows, suffices to power inte rstel- 
lar turbulence iVazguez-Semad eni. Pass ot. & PouauelllT99^ . The 
anticipated outer scale of the turbulence is around 100 pc (e.g., 
iBecket alJl993) . This is large compared with the size of individ- 
ual supernova explosion sites. Some supernova remnants can stay 
reasonably coherent for scales up to several tens of parsec. Fur- 
thermore, supernova explosions are known to be clustered, form- 



ing thereby large expanding patches called superbubbles (e.g., 
'Norman & Ikeuchi 1989: Ferriere 1992). Their scale is often larger 
than the scaleheight of the galactic disc, and they are primarily re- 
sponsible for driving fountain-like outflows. In any case, it appears 
that the localized point-like supernova explosions are able to act as 
a driver of the turbulence on a much larger scale. 

An obvious question concerns the connection between the ar- 
tificial forcings assumed in many computer simulations and any of 
the more realistic types of forcing. In the present paper we con- 
sider the properties of turbulence driven by random, localized ex- 
pansion waves, mimicking certain aspects of supernova-driven ex- 
plosions. One of the questions we are able to address with such a 
setup concerns information about the scale of the energy-carrying 
eddies of turbulence that is driven by small localized expansion 
waves. Is it true that such localized energy injections produce ki- 
netic energy predominantly at the scale of the original expansion 
wave, even though this scale is in general quite small? What are 
the effects, if any, that could shift the dominant power to larger 
scales. Large scales might plausibly arise if the duration of energy 
injection at one particular site becomes comparable to or larger 
than the turnover time of the resulting turbulence. This would re- 
inforce the original expansion wave until it has grown to a larger 
size. These issues are not directly connected with compressibility 
or with the effects of cooling or the nature of the equation of state 
of the gas. Thus, although interstellar turbulence is certainly highly 
supersonic, we consider here the case of subsonic turbulence in or- 
der to make the connection between solenoidal forcing in the form 
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of plane waves and irrotational forcing driving localized expansion 
waves. 



2 METHOD 

We solve the compressible Navier-Stokes equations assuming a 
polytropic equation of state relating the pressure p to the den- 
sity p via p — Kp'' with 7 = 5/3 and polytropic coefficient 
K = Csq/(7Pq^^), where Cso and po are constants. The momen- 
tum equation can then be written in the form 

Du _ 
where 



-V/l + / + fviBC, 



(1) 



h : 



PoCsO 
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is the enthalpy, D/Dt — d/dt + u ■ V is the advective derivative, 
and fvisc = P^^V ■ {2pvS) is the viscous force, where Sij — 
I J + — V ■ w is the traceless rate of strain tensor. The 
density obeys the continuity equation. 



Dlnp 
Dt 



-V ■ u. 



(3) 



We adopt a gaussian potential forcing function / of the form 
f{x,t)^V<t> (4) 
with 



,t) = N ei^p {[x - X{{t)f / R^] . 



(5) 



where x = {x,y,z) is the position vector, xt{t) is the random 
forcing position, R is the radius of the gaussian, and A'^ is a nor- 
malization factor. We consider two forms for the time dependence 
of Xf . First, we take Xi such that the forcing is 5-correlated in time. 
Second, we include a forcing time 5fforco that defines the inter- 
val during which Xf remains constant. On dimensional grounds the 
normalization is chosen to be A'^ = Cso \/ Cso-R/At, where 

At = max {5t, Stfovcc) (6) 

is the length of the time step, 5t, in the 5-correlated case or equal 
to the mean interval (Jtforcc during which the force remains un- 
changed, depending on which is longer. We begin by consider- 
ing the nature of flows generated in each case in the absence 
of a magnetic field. We use the PENCIL CODE,^ which is a 
non-conservative, high-order, finite-difference code (sixth order in 
space and third order in time) for solving the compressible hydro- 
dynamic and hydromagnetic equations. We adopt non-dimensional 
variables by measuring speed in units of the sound speed, c^o, and 
length in units of l/fci , where fci is the smallest wavenumber in the 
periodic domain. This implies that the nondimensional size of the 
domain is (27r)'^. 



3 RESULTS 

3.1 Delta-correlated forcing function 

We begin with the case where (5fforco ~ 0, so the forcing function 
is 5-correlated in time. In FigsQand|2|we show the resulting time 

|http : / /www. nordita ■ dk/software/pencil-codel 




Figure 1. Time averaged energy spectra for R = 0.2 (dotted line). 0.5 
(dashed line), and 1 (solid line). The durations of the runs are A'^tum = 333, 
224, and 81 turnover times, respectively. 128"^ mesh points, u = 10"'^, 
and <5-correlated forcing. The vertical arrows indicate the effective forcing 
wavenumber, as obtained from Eg. Isl. 



averaged energy spectra for two different resolutions. The results 
are relatively robust with respect to changes in the resolution. Also, 
the spectra are found to converge quite quickly after only about 5- 
7 turnover times, suggesting that there are no strong transients. For 
the results presented below, the durations of the runs, Tjun, is given 
in terms of the turnover time. 



urn — ^rms^pcakTrun 



(7) 



where fcpeak is the effective forcing wavenumber which, in turn, 
depends on the radius of the initial expansion wave, R. Indeed, the 
energy spectra show maximum power at the wavenumber fcpeak, 
which is approximately inversely proportional to R, with 



•-peak 



R : 



(8) 



The wavenumber fcpcak agrees with the width of the Fourier- 
transformed forcing function. [Note that the Fourier transform of 
is exp( — fc^i?'^/4).] This result is therefore not sur- 
prising, but it confirms quite clearly that a (5-correlated forcing 
function consisting of localized expansion waves produces maxi- 
mum energy injections at the scale of the expansion wave itself. 
This would therefore not explain the energy spectrum seen for ex- 
ample in simulations of supersonic turbulenc e (Porter et al. 199^ 
IPadoan & Nordlun(il99^lHaugen et all2004 or those inferred for 
interstellar turbulence faan et alj2004 . where energy injection is 
assumed to occur at the scale of the domain. 

The energy spectrum shows a dist inct subrange. Such a 
slope is predicted for shock turbulence iKadomtsev & Petviashvilil 
Il973h . although here the flow is subsonic and without shocks. 
However, it may be interesting to note that a spectrum has 
also be en seen in thejrr otational component of transonic turbu- 
lence jPorter et alJll998l) . Energetically, the irrotational compo- 
nent is subdominant compared with the solenoidal components 
(Porter et al. 1998; Padoan & Nordlund 1999; Hansen et al. 2QQ%. 
In the present simulations, however, there is hardly any vorticity 
and the flow is entirely dominated by the irrotational component. 
This, in turn, is a consequence of the irrotational nature of the forc- 
ing function combined with the fact that the viscosity is low, so that 
very little vorticity is produced. We consider vorticity production in 
more detail in Sect. l3.3l below. 
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Figure 2. Same as Fig.0 but for 256^ mesh points, f = 5 X 10 for Figure 4. Dependence on viscosity. Time averaged energy spectra for v = 

R = 0.2 (dotted line), 0.5 (dashed line), and 1 (solid hne). The durations 2 X lO""* (solid line), = 5 X IQ-"* (dashed Hne), u = lO"^ (dotted 

of the runs are A'^turn = 120, 64, and 29 turnover times, respectively. line). The durations of the runs are A'^turn = 28, 17, and 138 turnover times, 

respectively, and the values of the St are 0.23, 0.19, and 0.17, respectively. 

■Jiforcc = 1, B. = 0.2, 128^ mesh points. 




Figure 3. Time averaged energy spectra for R = 0.2 for 5tforco = 1 
(St = 0.08, dashed line), ^tforcc = 0.2 (St = 0.05, solid line), and 
^itoicc = (i.e. (5-correlated, St = 0, dotted line). The durations of the 
runs are Nt^vn = 16, 64, and 333 turnover times, respectively, and the 
rms velocities are 0.008, 0.027, and 0.058, respectively. 128^ mesh points. 



3.2 Forcing function with memory 

Next we look at the case where the expansion wave persists for a 
certain amount of time, (5fforco- In Fig.|3|we show time averaged 
spectra for different values of (5f force- This time interval is conve- 
niently expressed in non - dimensional form as a Strouhal number 
Landau & Lifshitzll987tlKrause & RadleJl98(j) . 



St = Sttoi 



(9) 



Note that for larger values of Sttovcc the spectrum consists of sep- 
arated bumps, but the overall power appears to be shifted to larger 
scales and smaller wavenumbers. The presence of bumps is sug- 
gestive of an acoustic resonance or standing wave pattern that may 
be excited in the simulation box. Such a phenomenon is unfamil- 
iar in the context of mostly vortical turbulence. It is plausible, and 
indeed compatible with our as yet limited results, that these bumps 
can only occur if there is sufficient scale separation between fcpcak 
and ki. Furthermore, bumps in the spectrum seem to be possible 
only when the memory time is not too short, i.e. when the effects 
of temporal randomness are limited. 




Figure 5. Dependence on resolution. Time averaged energy spectra for 
R = 0.2, using 512^ mesh points (solid line, u = 5 X W'''), 256^ 
mesh points (dashed line, u = 10~*), and 128"^ mesh points (dotted line, 
1/ = 2 X 10—^). The durations of the runs are A^turn = 22, 59, and 28 
turnover times, respectively, and the values of St are 0.25, 0.24, and 0.23, 
respectively. The insets shows the evolution of Mrms and Wrms for the run 
with 512"^ mesh points. 



In Figs|4]and|5]we demonstrate that decreasing viscosity and 
increasing resolution produces more power at small scales, gradu- 
ally building up a fc^^ power law between the forcing scale and the 
dissipative cutoff scale when the resolution is sufficient for those 
scales to be resolved. 

Visualizations of the density (Fig. |6} show how the initially 
highly ordered expansion waves turn rapidly into a complicated 
pattern.^ By the time t — 50 the velocity has reached a statisti- 
cally steady state and the flow has ceased to bear any resemblance 
to the initial expansion wave. 



^ Animations of the density can be found at 
|http : / /www ■ nordita ■ dk/ "brandenb/movies /gauss ianpot| 
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t= 2^0 



t= 5^0 



'(3) 





t= 100^0 




Figure 6. Visualization of In p on the peripheiy of the box at different times, for kiR = 0.2,Re = 50, St = 0.25, u = 5 X 10 ^, and 512^ mesh points. 
Note that in the fully developed state individual expansion waves can hardly be recognized. 



3.3 Vorticity production 

One of the main differences between the present forcing function 
and those used in some previous p apers ( Eswaran & Pope 1988; 
IPorteretal. 1998: Padoa n & Nordlund 2002: Haugen et al. 2004) 
is that here the forcing is completely irrotational, i.e. V x / = 
and V ■ / 7^ 0. Therefore, unlike previous cases where instead 
V X / 7^ and V • / = 0, and where vorticity (u> = 'V x u) was 
immediately produced, here vorticity can only be produced through 
viscous interactions. This can be seen by writing the viscous force 
in the form 



^ VlSi 



~uV X V X u + |!^VV 



(10) 



where d = SijVj In p is a term that could drive vorticity even 
with zero vorticity initially. Indeed, the vorticity equation, obtained 
by taking the curl of Eq. Q, can then be written as 

'dt 



= V X (u X u;) + vV ijj + vV X G. 



(11) 



Note the absence of the baroclinic term, i.e. the cross product of 
temperature and specific entropy gradients. This is because we have 
restricted ourselves to a gas with an isothermal equation of state. 
The b aroclinic term is known to produce vorticity (e.g. Korpi et al. 
Il999b ). The nonlinear term, V x (u x tj), could in principle lead 
to exponential amplification of an initial seed vorticity — much like 
a term in the small scale dynamo problem, discussed in the next 
section. However, no evidence for such an effect has been found in 
the present simulations. Thus, the only remaining term is V x G 
which is due to viscosity. One might therefore expect that as we 
decrease the Reynolds number, defined here in terms of fcpcak as 



Re 



(12) 



the amount of vorticity production decreases. The numerical simu- 
lations yield values for the normalized vorticity, 



normalized vorticity — ujrms/ 



peak } 



(13) 



that are essentially within the 'noise'. Indeed, the data shown in 
Table □ are well below unity and do not give a clear trend as a 
function of resolution or Reynolds number. For all these runs the 
Strouhal number is St ~ 0.25. 

In conclusion, the amount of vorticity production depends 
mostly on the numerical resolution, suggesting that no measurable 
vorticity is produced by physical effects. 

3.4 Dynamo action? 

In view of the formal analo gy between the vorticity and induction 
equations jBatchelodfTisol) . it appears that dynamo action should 



Table 1. Results for the normalized vorticity, t<Jrms/(Mrmsfcpcak)i as a 
function of Re and resolution. The durations of the runs vary between 



16 and 250 turnover times, and u varies between 10 



and 



5 X 10 ^. For the high resolution run with 512^ meshpoints and u 
5 X 10~^ we have A'^turn = 59 turnover times. 



Re 



512^ 



256=^ 



128^ 



50 1.4 X 10-3 



25 
12 
4 



8.7 X 10-3 
1.1 X 10-2 
1.6 X 10-2 
7.6 X 10-3 



2.9 X 10-3 
2.0 X 10-3 
1.5 X 10-3 



be as difficult to achieve as the amplification of vorticity by the 
nonlinear term, V x (u x cj). In order to address this problem, we 
solve the induction equation for the magnetic field B, 



V X (u X B) +r]V B, 



(14) 



dB 
'dt 

simultaneously with the other equations. Equation Q gains a 
term — the Lorentz force per unit mass, J x B / p, but this effect 
would only become important once the magnetic field becomes 
strong. (Here, J = V x B / fio is the current density and fj,o is 
the vacuum permeability.) As initial condition we adopt a magnetic 
vector potential that is 5-correlated in space. This results in a mag- 
netic energy spectrum proportional to fc'*. 

In all the cases that we have investigated so far, we have 
not found evidence of field amplification, i.e. there is no dynamo 
action. Our highest resolution run has 5123 jjjggjj points and a 
magnetic Reynolds number, Rm = Wrms/(??fcpoak), of 250 and 
Re = 50. Kinetic and magnetic energy spectra are shown in Fig.0 
For this run the magnetic decay rate, A = —din Bmis/dt, is about 
30 times the ohmic decay rate, ?7fcpcak- This is approximately equal 
to the scale separation ratio between fcpcak and the dissipative cut- 
off wavenumber. Such a rapid decay rate suggests that this type of 
flow has a highly destructive effect on the magnetic field. Energy 
spectra of the magnetic field confirm that the spectral energy decays 
at small scales, while it stays unchanged at large scales. 

The absence of any dynamo action agrees with earlier re- 
sults for purely irrotational turbulence liBrooka»1999i) . On the other 
hand, earlier analytic considerations suggested that dynamo action 
should be possible for irrotational flows and that the growth rate 
should increase with incre asing Mach numb er to the fourth power 
(Kazantsev et al. 1985: Mo ss & Shukiirovll9 96V Subsequent work 
by H augen et aL ^200 4 ) did not discuss growth rates, but critical 
magnetic Reynolds numbers which shows that the small scale (non- 
helical) dynamo becomes about twice as hard to excite when the 



Random localized expansion wave turbulence 5 




Figure 7. Kinetic and magnetic energy spectra, E(k) and M(k), re- 
spectively, at three different times, separated by At = 205tfoi.ce> with 
for R = 0.2, Rni = Mrms/(r?fcpeak) = 250, Re = 50, and 512^ 
mesh points. The duration of the run is A'turn = 22 turnover times. 
iy = 5 X 10-5,17 = 5 X 10-5. 

Mach number exceeds unity. This was explained by the presence of 
an additional irrotational contribution in the supersonic case. This 
contribution remained however subdominant and did not contribute 
to the dynamo. This is consistent with the present result that up to 
the largest magnetic Reynolds numbers accessible today, purely ir- 
rotational turbulence does not produce dynamo action. 



4 CONCLUSION 

The localized expansion waves used as forcing functions in the 
present paper resemble in some ways the driving by supernova ex- 
plosions. An obvious difference is that the supernova driving pro- 
duces highly supersonic flows. This may be important in connec- 
tion with understanding the amount of vorticity production in irro- 
tationally forced flows. Indeed, Koroi et al. ( 1999b) found that the 
main driver of vorticity production in their supernova driven turbu- 
lence simulations is the baroclinic term. 

The present work suggests that the turbulence would indeed be 
peaked primarily at small scales. Memory effects (i.e. finite values 
of (5f force) contribute somewhat to enhancing power at larger scales, 
but, because of the dip in power at intermediate scales, the resulting 
energy spectrum is still very different from the power law spectra 
obtained with forcing directly at large scales. The main difference 
is the absence of vorticity in the present simulations. Since vorticity 
is mainly driven by the baroclinic term, we may expect that it is also 
this term that would primarily contribute to enhanced power at low 
wave numbers, corresponding to scales larger than the radius of the 
original expansion waves. 
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